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1. Introduction 

Bessel function values I n (z),J n (z) generated by BESLCI [1 J were compared with check values 
calculated via the ascending series 
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using multiprecision arithmetic [2]. 2 The 60-bit mantissa of each part of the test value was sub- 
tracted from the corresponding mantissa of the check value, and the result expressed as a multiple, 
ra, say, of the last bit. For most errors mentioned in this certification, | m | ^ 7. Bit comparison was 
used in preference to calculation of relative error in order to simplify computations. An error of 
m bits in the mantissa corresponds to a relative error between m • 2~ H() and m • 2~ 59 . 

This test is too strict, however, near a zero of the real or imaginary part of the function being 
tested. For complex values of z, it is more realistic to compare the bit error with the greater of the 
real and imaginary parts, since the relative error in a complex number £+ irj is (S£ + iSr))!^-^-^ | , 
5f and 8r) denoting the respective errors in the real and imaginary parts. For arguments in the 
neighborhoods of zeroes of Re Jn(z), Im J n (z), Re /«(z), or Im I n (z), an assessment of absolute 
errors was made to 60 binary places by right-shifting the mantissa of both test and check values 
before subtraction. This was done whenever the order n was less than \z\ and the modulus of the 

check value less than -. These cases are distinguished by asterisks in the computer printout; 

statistics for this test are compiled separately. 

A bit comparison test was also used to check the section of code involving the two-term as- 
cending series for small \z\. 

A second feature of BESLCI which was tested systematically is the error return. This portion 
of the code determines whether all orders /i=0, 1, . . ., NB~1 can be calculated to the desired 
accuracy. A modification occurs if NB is so large compared with \z\ that there is underflow in the 
representation of J m -\{z) or / N B-it). When this happens, the program determines the largest 
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value NCALC (> |z|) such that accurate values can be calculated for n = 0, 1, . . . , NCALC-1; 
these values are then calculated. To test this, BESLCI was called with NB=5600. In each case, 
the program returned NCALC < NB, insuring that the error return loop had been executed. Then 
the values of NCALC and the Bessel function of order NCALC-1 were checked by calling BESLCI 
with NB set equal to NCALC; this involved another execution of the error return loop. In each case, 
the new NCALC was identical with the old. Finally, the Bessel function value of order [|z|] was 
checked by calling BESLCI with NB = [ | z | ] + 1; this involved no execution of the error return loop. 
In both cases, bit comparisons on 60 significant bits were performed. 
Results of the tests are described in the next two sections. 

2. The Bit Comparison Test 

Arguments z were chosen as follows. Let Rj denote the set of rational numbers 2 J -p, where p is 
a 60-bit fraction, with \ ^ p < 1. Two complex arguments were generated in each of the sets 

Aj = {x + iy: xeRj , y = 0} j =-6(1)6 

B k ={x + iy:x = 0, yeR k } k =-6(1)6 

C jk .= {x+iy:xeRj, yeR k } ;, A --6(1)6, 

using a pseudorandom number generator [3]. One of the two arguments was used to generate /'s 
and the other J's. For each argument, functions of order 0, 1, . . .,10 were tested. It should be 
noted that the sets Cjk include complex arguments with real and imaginary parts of quite disparate 
magnitudes. 

Summary statistics for the tests are: 

Test on 60 Significant Bits 

903 results correct to 60 bits* 

1329 results in error by 1 multiple of last bit 

967 results in error by 2 multiples of last bit 

657 results in error by 3 multiples of last bit 

404 results in error by 4 multiples of last bit 

333 results in error by 5 multiples of last bit 

247 results in error by 6 multiples of last bit 

208 results in error by 7 multiples of last bit 

1422 results in error by more than 7 multiples of last bit 



6470 (Total) 

Test on 60 Binary Places 

320 results correct to 60 places** 

479 results in error by 1 multiple of last place 

279 results in error by 2 multiples of last place 

160 results in error by 3 multiples of last place 

82 results in error by 4 multiples of last place 

59 results in error by 5 multiples of last place 

45 results in error by 6 multiples of last place 



*Excluding 202 results for pure real or imaginary argument, for which test and check values are both zero. 
**Excluding 370 results for pure real or imaginary argument. 
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26 results in error by 7 multiples of last place 

88 results in error by more than 7 multiples of last place 



1538 (Total) 

In the first test, the maximum error was (6651)h. This occurred in the real part of Jn(z) at 
n=10 and z= (.1216 . . .) -£(.01928 . . .), for which /„(z) = -2- 71 -p-i • 2" 62 • q, with 
\ ^ p, q < 1. It is appropriate to assess this error relative to \Jn(z) \ . This can be done fairly simply 
by dividing the bit error by 2 9 , since 9 is the difference between the binary exponents of Re J ,,(z) 
and Im J n (z). This "relative error" is 7. On applying this technique to all the bit errors, the maxi- 
mum relative error was found to be (34)k. Allowing for the fact that one mantissa can be no more 
than twice the other, we see that the relative error is bounded by (70)s, showing that BESLCI 
yielded values correct to 54 bits out of 60. 

In the test on 60 binary places, the maximum error was found to be (77)s, showing that at 
least 54 binary places are correct. (Indeed, the first 57 places were correct in 94% of the cases.) 

The test on the two-term ascending series used only the argument set Cjk (see above), where 
j and k each take the values — 16 and — 100. The relative error was bounded by (60) 8 , in comparison 
with (70)8 for the test described above. 

3. Testing of Error Return 

The following arguments z = x + iy were used: 

V+i-2* j, A=-9(3)12 

V j =-9(3)12 

i-2 k k =-9(3)12; 

the values of NB are described in section 1. After them's were calculated using x+iy as the argu- 
ment, the /'s were calculated using y-\- ix. This showed that BESLCI gave results in accord (to 
60 bits) with the equations 



Jn(z)=Jn(z), In(z)=(-i) n J n (iz). 

In the test on Bessel function values of order NCALC-1, the largest error was 5 multiples of 
the last bit; most were zero. 

In the test on order [ | z | ] , the errors were larger, owing to the greater number of back-recursion 
steps. To calculate ./jow (4096 + i • 2~ 9 ), for example, the subroutine call with NB = 5600 involved 
approximately 900 more steps than the call with NB = 4097. Bearing in mind that each step includes 
2 multiplications and 3 additions for both the real and imaginary parts of the function value and 
also that the computer uses truncation in floating-point arithmetic rather than roundoff, we see 
that the error of (2234) H = 1 180 in the last bit in the real part of/ 4 096 (4096 + i • 2~ 9 ) is to be expected. 
The imaginary part of J wm (4096+ i ■ 2~ 9 ) was about 2~ 15 times the size of the real part, and the 
bit error, (72461407) H , was about 2 14 times as big, which, again, is to be expected. 

The errors mentioned in the preceding paragraphs were the worst cases. For all other argu- 
ments, the errors were within the expected tolerance, based on the number of recursion steps, and 
the relative sizes of x and y. 

4. Summary and Conclusions 

For /7 = 0(1)10 and the values of z tested, BESLCI yielded J n (z) and I„(z) on the UNI VAC 
1108 with either a relative error or an absolute error bounded by 10 -16 . The latter assessment applied 
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when n < \z\ and the magnitude of the real or imaginary part of the Bessel function value being 
tested is less than i; the former applies in all other cases. 

The error return feature was tested and found to be executed accurately in all cases. This test 
also confirmed the error bounds described above for values of z in the square | Re z\ < 64 and 
| Im z | < 64. 

For other values of n and z, similar accuracy may be expected since the tests we have used 
include large and small values of \z\ and \z\/n. 

Calculations were performed on a UNI VAC 1108 under EXEC 2. All were carried out with 
double-precision mantissa of 60 bits, corresponding to just over 18 figures. Subroutine executions 
took approximately 0.5N milliseconds, where 

7V=max (|z|, NB) + 10. 
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